Classical Nuclear Motion in Quantum Transport 
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An ab initio quantum-classical mixed scheme for the time evolution of electrode-device-electrode 
systems is introduced to study nuclear dynamics in quantum transport. Two model systems are dis- 
cussed to illustrate the method. Our results provide the first example of current-induced molecular 
desorption as obtained from a full time-dependent approach and suggest the use of ac biases as a 
way to tailor electromigration. They also show the importance of non-adiabatic effects for ultrafast 
phenomena in nanodevices. 
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The study of electron-nuclei interaction (ENI) in solids 
has a long and important history, dating back to the early 
days of quantum mechanics. Since then, the ENI has be- 
come a pillar concept of our understanding in condensed 
matter [l|. Increase in computer power and the intro- 
duction of ab-initio molecular dynamics |2| have made 
possible quantitative studies of ENI in many materials. 

Recent advances in nanotechnology pose new ques- 
tions about the ENI in out-of-equilibrium open systems 
at the nanoscale Among the techniques used to 

study the ENI in this context, quantum transport ex- 
periments stand out as a special case, since charge con- 
duction is at the same time a way to characterize the 
nanodevice and a property to be exploited in its oper- 
ating regime. A case in point is molecular junctions, 
where electron injection can stimulate local vibrations 
0, and possibly electromigration Q . More generally, 
describing ENI in time-dependent (TD) quantum trans- 
port is expected to become of great technological interest, 
since future nanodevices will operate under the influence 
of ever faster time varying external fields. Accordingly, 
those regarded at present as marginal transient effects 
will soon become center stage features. Assessing and 
engineering the ENI is then a key ingredient to increase 
the device efficiency 111- To date, most theoretical 
studies have addressed the ENI in steady-state phenom- 
ena H [H El E E H El and often perturbatively in 
the nuclear displacements Hfij . Going beyond the 

harmonic approximation and including the ENI in a first- 
principles, TD framework is a difficult theoretical task 
which has received so far scarce attention [lH \vk Hil [Tq| . 

In this work we propose a first principle approach to 
TD quantum transport which treats the nuclear mo- 
tion in the Ehrenfest dynamics (ED), and the elec- 
trons within Time-Dependent Density Functional The- 
ory (TDDFT)[lj. The ED has been extensively used 
in several contexts; in quantum transport, it was con- 
sidered in [l^. The ED correctly displays many impor- 
tant features of ENI, but gives an incomplete account 
of the Joule heating of the nuclei by the electrons [lfSj . 
However, such effect is not the aim of this first work. 
Here we use ED which, while treating ENI at the mean 



field level, includes, as opposed to the Born Oppenheimer 
Approximation, electronic transitions. Our approach to 
transport permits the description of transient effects: it 
can be used to describe history dependent currents, hys- 
teresis phenomena, etc. Another main advantage is that 
it can perform the TDDFT-ED of devices connected to 
infinitely long leads. 

We illustrate the approach with two model devices, 
where the electrons interact only via the ENI: a Holstcin 
wire (Dl) and a diatomic molecule (D2). We choose these 
two rather different systems to show the versatility of our 
method in discussing transient phenomena and overcom- 
ing the limitations of adiabatic treatments. More in de- 
tail, our results show that: (i) For weak electron-phonon 
(e-ph) couplings, Dl exhibits an almost periodic nuclear 
displacement (with period=l/density), reminiscent of a 
Peierls distortion. On applying a dc bias, the nuclei oscil- 
late (with decreasing amplitude) and the period changes 
to accommodate the current flow. On increasing the e- 
ph coupling, Dl changes from conducting to insulating, 
(ii) D2 is deformed by a small, suddenly switched on dc 
bias. Above a critical value of the bias, the molecule dis- 
sociates. This is the single most important result of this 
work. To our knowledge, it provides the first example of 
current-induced molecular desorption as emerging from 
the full TD dynamics of a nanodevice. (iii) The des- 
orption cannot be described in any adiabatic formalism, 
since it is due to electronic excitations induced by the 
nuclear motion, (iv) The desorption can be tuned by the 
intensity and frequency of an ac bias, suggesting a way 
to control electromigration in molecular devices. 
The Method. Following Refs. [U |H Hj|, in the ini- 
tial ground state the central region (C) is in contact to 
seminfinite left (L) and right (R) leads. With ENI, the 
Hamiltonian for C reads 



M 

ff c [x]= v ^Mc 3 , 



(1) 



where M is the number of one-electron states of C and 
x = (.Ti, . . . , xn) are the nuclear coordinates. Outside C, 
the nuclei are clamped. In the Kohn-Sham (KS) scheme 
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of TDDFT 20 . Vij would be the matrix element of 
the KS Hamiltonian. The nuclear classical Hamiltonian 



is H cl = Y,k=iPk/( 2rn k) 
system is governed by 



C/ C i(x). The dynamics of the 
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where |v?(i)) is the many-electron state at time t, and 
i? c i[x] is the electron Hamiltonian of the contacted sys- 
tem L+C+R. Given a configuration x, flei[x] is a free- 
particle Hamiltonian. The many-electron ground state 
|^ 9 [x]) consists of bound, resonant, fully reflected waves, 
plus left and right going scattering states. The paramet- 
ric dependence of H e \ on x renders every eigenstate a 
function of x. The ground state value x = x s is com- 
puted using a damped ground-state dynamics: starting 
from an initial x the coordinates are evolved according 
to m k x k = ~jx k - d [U d + (* fl |i?c|* ff )) /dx k , with 7 
the friction coefficient. 

Having x g and the corresponding {-ip} one-electron or- 
bitals, we apply an external bias and evolve the system. 
Assuming metallic electrodes and instantaneous screen- 
ing, the size of C is chosen so that the potential drop 
at any time occurs entirely in C. Thus, the TD part 
of the electrode Hamiltonian is a spatially uniform shift 
Ur)(t),r] =L,R. We use a novel mixed quantum-classical 
evolution algorithm, which combines a recently proposed 
generalization of the Crank-Nicholson method 23] for the 
{ijj} with a Verlet-like integrator for the x. Schematically, 
in terms of the discretized time t m = 2mS, 



|^(m+l)| = {S{t m , X <- m )}lp( m )} 



(toH 
Pk 



(m) 
= Pk 



I Pk 



(m+ 
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+ 25F k [x < - m \{iP (m +^}} 



45p k m+1) /m k 
+ 25F k [x 
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(m+2) 



with ^("0 = ^(t m ), x (™) = x(f m ) and p( m ) = p(* m ). 
The unitary matrix S depends on time through x(£) and 
the TD bias U v (t), rj =L,R. Full details of the electronic 
evolution can be found in Ref. |2^|. The force F[x, {ip}} 
is given by the r.h.s. of Eq.©. 

To illustrate the method, we describe electrodes L and 
R in terms of one-dimensional tight-binding (TB) Hamil- 
tonians with a hopping V between nearest neighbor sites 
(TB parametrization of Au wires suggest |V| in the range 
0.3 -7- 1.0 eV). Left and right going scattering states have 
energy e within the band (— 2|V|,2|V|) and can be ob- 
tained by solving the Schrodinger equation in C with ap- 
propriate boundary conditions. Bound state eigenener- 
gies Eb < — 2\V\ satisfy Det [e&l — Hq — E(eji)] = 0, and 



the associated wavefunction is given in C by the kernel 
of [e&l — Hq — (ffc is the projection of Hq onto 

the one-electron Hilbert space, is the embedding 

self-energy). The topology of region C might also lead 
to states rigorously confined in C (see below). These are 
resonant eigenstates of the uncontacted He with zero 
amplitude at the interface with the two electrodes. 
Model device Dl. The semiclassical Holstein model is 
a valuable tool to gain insight into e-ph interactions 
[iil Ei^j ] . Here, we investigate TD transport through a 
Holstein wire described by 



H C = V 



M-l 

E 

i=-M 



M 



c\c i+ i + h.c. - g 



^ Xifii, (5) 

i=-M 



where hi — c\ci is the local density operator and Xi are 
the phonon coordinates. The nuclear classical potential 
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FIG. 1: Left panel: Ground state displacement Xi, i — 
— 7, at weak (A = 0.1) and strong (A = 1.0) coupling. 

The inset shows the LDOS (energy in units of |V|) on site —6 
(dashed line) and site (solid line). Right panel: current I(t) 
(in units of 10 -2 |V|) at weak and strong coupling along the 
bond (-7,-6) (solid thin), (0,1) (solid thick) and (6, 7) (dashed) 
after the sudden switching on of a bias Ul = 0.5|V| in elec- 
trode L. The insets show a 3D plot of x(i) between t = and 
t — 480|V| _1 (in the top panel the range of x is between -3.2 
and -4 while in the bottom panel is between -19.4 and -20). 



is U c \(x) = ^Ei=-m ^o 1 !- The strength of the e-ph 
interaction is determined by the dimensionless parameter 
A = g 2 /(2Vuio)- In the original derivation by Holstein, 
Xi corresponds to an internal coordinate (bond-length) at 
the i-th site. We study Dl at half-filling (ep = 0), in the 
adiabatic regime (a — ujo/V = 0.1) at weak (A = 0.1) 
and strong (A = 1) coupling. A damped ground-state 
dynamics, as described above, was used to get the con- 
verged ground-state (Fig. ^ left) f° r a region C with 
M = 7. The energy spectrum between — 2\V\ and £f 
was discretized and good convergence was achieved with 
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Nk = 500 meshpoints. For A = 0.1, a Peierls-like distor- 
tion is seen: an even-odd behavior of x\ is manifest (in 
general, P — l/n, and here n = 0.5), but exact periodic- 
ity is prevented by the finite size of C. The inset shows 
the local density of states (LDOS) close to the interface 
and in the center, respectively. At e < —2\V\ we ob- 
serve three peaks due to three bound states. The picture 
changes dramatically at A = 1: The number of bound 
states equals the dimension of C and, almost uniformly, 
Xi ~ —20. Only the x% close to the interfaces are slightly 
above this value. 

In Fig. Q right, we plot the current lit) in three 
different points of C after suddenly switching on a bias 
Ul = 0.5|V^| in the left electrode. All calculations use a 
time step 6 = 0.01|V| -1 . For both A = 0.1, 1.0, the cur- 
rent in the (short) transient is similar to the case with the 
Xi clamped (not shown), since electrons are much faster 
than nuclei. At A = 0.1 the steady current shows super- 
imposed oscillations of frequency ujq. Instead, for A = 1, 
all the 2M + 1 = 15 bound states in C are occupied 
(including the Hartree potential would have led to a sig- 
nificant reduction of this excess density in C) , no current 
fluctuations occur in the center and at the right interface 
and only a very short transient is observed at the left in- 
terface. The insets display the TD Peierls distortion. For 
A = 0.1 all the x's oscillate with an exponentially decreas- 
ing amplitude (within our simulation time t = 480 1 | 1 ) . 
The overall shape of the Xi changes to accommodate the 
net electron flow. For A = 1 no current flows and only the 
Xi close to the left interface oscillate. The decay of the 
amplitudes of the x% is partially due to the inefficiency of 
ED in transferring energy from electrons to nuclei. We 
also considered the sudden removal of a bound electron, 
as obtainable for example by optical means: in this case 
Dl provides a strong transient oscillating response (not 
shown). On speculative grounds, such behavior could be 
used for ENI based photosensors. 

Model device D2. We consider a central region C with 
the simplest non-trivial topology, a four-atom ring. This 
is our model molecular device D2, with nuclear positions 
= (x i: yi), i = 1, .., 4 (Fig. |3 top-left). In D2, only nu- 
clei 2 and 3 (N2,3) are let to move in the xy plane. The 
origin of the xy plane is the midpoint of nuclei 1 and 
4. For the hopping parameter we choose the form j2|| 
Vi^j = V c ^^-, Vi=j=0, where r y - = |r; The purely 

repulsive classical term is given by U c i = \ 2i=y ^/ r ij- 
We tune the parameters V c , A and ru to ensure a reason- 
able domain of structural stability for D2. For V c = 4V, 
A = 0.75|V^| and ru = 2, D2 is stable against deforma- 
tions up to ~ 10% of the equilibrium distances. Also, 
we get the ground state positions r 2 = (0,0.8313) and 
r*3 = — V2 (the 2 — 3 symmetry remains true in the pres- 
ence of the bias). The LDOS has one peak below — 2|V|, 
i.e. one bound state. There is also a resonant state 
|^) = -4=(|2) - |3)) with energy \V 2 z\ > inside the band. 
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FIG. 2: Ih,n(t) (in units of \V |) at the left (solid line) and right 
(dashed line) interface for a sudden switching Ul = 0. 25 1 V 
(green[light grey]), 0.5|V| (red [thin dark grey]) and 1.0[V| 
(blue [thick black]). Inset a): Time dependent density of 
nucleus 3 (by symmetry ri2 = nz at any time). Inset b): xz(t) 
(solid line) and j/3 (t) (dashed line) of nucleus 3 (by symmetry 
X2 = 3:3 and 2/2 = — 2/3 at any time). Time scale and color 
coding in insets a) and b) is the same as in the main hgure. A 
schematic of the device D2 and the adiabatic result (labeled 
curve) are also shown. 



The system is taken at half filling, and we switch on a 
bias U-l at t = 0. As for Dl, 6 = 0.01|V| _ \ N k = 500; 
for the masses of N2,3 we choose lOOjX^I" 1 . 

In Fig. [21 we plot the TD currents 7l and Ir at 
the left and right interfaces, respectively. At small bias 
U L = 0.25|V|, Il.r rapidly increase and after a short 
transient (with the nuclei essentially still) start to os- 
cillate around a steady value. Inset b) shows the cor- 
responding nuclear dynamics. The equilibrium rhom- 
bic geometry changes and the molecule gets deformed 
in the biased system, with N2,3 having damped oscil- 
lations around two new positions (for the damping and 
ED, considerations similar to the case of Dl apply). We 
also notice from inset a) that the charge density of N2,3 
slightly increases. 

Highly interesting is the strong bias case U L = 1.0|V|. 
After a very short transient, Il,r sharply decrease (rather 
than oscillating around a steady value) and become zero 
at to ~ lOOl^! -1 . After to, 7l,r separate: 2l increases 
while Jr decreases, reaches a negative minimum and then 
increases to eventually re-join Zl at t\ ~ 160| —1 . For 
t > t\, Jl — Ik and their value equal the steady cur- 
rent (calculated from the Landaucr formula) of the chain 
without N2,3. The behavior of Il,r can be understood 
looking at the nuclear dynamics (inset b): the force ex- 
erted by the electron flow is strong enough for atomic 
migration to occur. N2,3 are pushed to the right by the 
current, overcome the confining potential and get disso- 
ciated from region C. Thereafter, they form a diatomic 
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FIG. 3: Left panel: / L , R ({) (in units of |V|) at the left (solid 
line) and right (dashed line) interface for Ul = 1.0|V|,w = 
0.05|V|. In the inset, X2(t) = xz(t) (thick solid) and yz(t) 
(thin solid), yz(t) (dashed). The time scale is the same as 
in the main panel. Right panel: same as left, but with lu — 
0.10|V|. Left and right insets have the same vertical scale . 



molecule vibrating along y (see inset b) and traveling 
along x at uniform speed. The pronounced minimum of 
Jr shortly after t is due to a sudden charge transfer from 
electrode R (via atom 4) to the diatomic molecule when 
N2,3 pass above nucleus 4. This is confirmed in inset a) 
where the density suddenly increases in correspondence 
of the minimum in Jr. We also note that the total den- 
sity of the dissociated molecule is about 1 (exact charge 
quantization would occur in the adiabatic approximation 
only). Here, including the Hartree potential would not 
lead to qualitative changes of the dissociation process. 
A more realistic model should also include the possibil- 
ity for the molecule to chemisorb onto the electrode, to 
fragment etc., but this is beyond the scope of the present 
Letter. 

To address the dissociation mechanism, we study two 
different Ui,(i), with the same asymptotic value 0.5| V| . 
For a sudden switching, one observes a behavior delayed 
but similar to the 1.0|V| case above. Instead, for an adi- 
abatic switching j^l the dimer does not dissociate ( in 
Fig. |21 the current is displayed in orange and labeled 
"adiabatic" ). We conclude that electronic excitations 
induced by the nuclear motion play a crucial role in the 
electromigration, a role that can not be accounted for 
by any adiabatic formalism (without transients the nu- 
clei experience no force). However, we observe that ED 
might overestimate the value of the critical bias for which 
dissociation occurs. 

We next examined the response of D2 when subject 
to a high amplitude ac bias Ut(t) = Ul sin(wt), with 
U L = 1-0|V|. At low uj (Fig. [3 left), the system qualita- 
tively behaves as in the dc clear from both cur- 
rent (main left) and coordinates (inset left) panels. The 
nuclei overtake the barrier before the change in Ul(£) 
produces a force able to "recall" them back. After the 
desorption, Jl,R oscillate as they would for a linear chain 
without N2,3. At larger to (Fig. |3 right) the dissocia- 
tion is delayed, since atoms 2 and 3 are recalled back a 



few times before leaving region C. Eventually, for high 
enough frequency (> 0.3 |V|, not shown) atoms 2 and 3 
oscillate without leaving C (within our simulation time) . 
These results point to a possible use of ac biases as a way 
to tailor molecular desorption in nanodevices. 

In conclusion, we presented a mixed quantum- 
classical scheme to describe electron-nuclei interactions 
in quantum transport. The scheme is straightforwardly 
amenable to an ab-initio implementation. Our results 
show the necessity to go beyond adiabatic schemes and 
to use a full time-dependent, nuclear-dynamics approach 
for ultrafast transient phenomena which are expected to 
become important in future generation nanodevices. We 
acknowledge useful discussions with G. J. Ackland, P. 
Bokes and E. N. Economou. This work was supported 
by the European Community 6th framework Network of 
Excellence NANOQUANTA (NMP4-CT-2004-500198). 
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